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Abstract. We discuss open problems related to the stochastic modeling of cardiac 
function. The work is based on an experimental investigation of the dynamics of 
heart rate variability (HRV) in the absence of respiratory perturbations. We consider 
first the cardiac control system on short time scales via an analysis of HRV within the 
framework of a random walk approach. Our experiments show that HRV on timescales 
of less than a minute takes the form of free diffusion, close to Brownian motion, which 
can be described as a non-stationary process with stationary increments. Secondly, 
we consider the inverse problem of modeling the state of the control system so as to 
reproduce the experimentally observed HRV statistics of. We discuss some simple toy 
models and identify open problems for the modelling of heart dynamics. 
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1. Introduction 

The human heart does not beat at a constant rate, even for a subject in repose. Rather, 
there is strong variability of the heart rate. The complexity of this heart rate variability 
(HRV) presents a major challenge that has attracted continuing attention. Many of 
the explanations proposed are by analogy with paradigms used in physics to describe 
complexity, including: deterministic chaos [I]; the statistical theory of turbulence [2]; 
fractal Brownian motion [3J; and critical phenomenon pE]. They have led to new 
approaches and time-series analysis techniques including a variety of entropies EJ [7] , 
dimensional analysis [8] , the correlation of local energy fluctuations on different scales 
[9], the analysis of long range correlation [10], spectral scaling [Til H2], the multiscale 
time asymmetry index [13J, multifractal cascades [TH [To] . All these measures allow one 
to describe HRV as a non-stationary, irregular, complex fluctuating process. Depending 
on the technique in use there has been a very wide range of conclusions about the 
regulatory mechanism of heart rate, ranging from a stochastic feedback configuration 
[T6] to the physical system being in a critical state [9]. HRV can also be considered in 
terms of the interactions between coupled oscillators of widely differing frequencies j!7j . 

Although we now have this huge variety of tools and approaches for the analysis 
of HRV, only the last-mentioned has enabled us to understand the origins of some of 
the time-scales embedded in HRV. Each time-scale (frequency) in the coupled oscillator 
model [T7] is represented by a separate self-oscillator that interacts with the others, 
and each of the oscillators represents a particular physiological function. The frequency 
variations in HRV can therefore be attributed to the effects of respiration (~0.25Hz), 
and myogenic (~0.1Hz), neurogenic (~0.03Hz) and endothelial (~0.01Hz) activity. 
HRV also contains a fast (short time-scale) noisy component which forms a noise 
background in the HRV spectrum and can be modelled as a white noise source [17] . 
Its properties are currently an open question, and one that is important for both 
understanding and modelling HRV. A practical difficulty in experimental investigations 
is the presence of a strong perturbation, respiration, that occurs continuously and exerts 
a particularly strong influence in modulating the heart rate. This modulation involves 
several mechanisms: via mechanical movements of the chest, chemo-reflex, and couplings 
to neuronal control centres [18]. Spontaneous respiration gives rise to a complex non- 
periodic signal, and this complexity is inevitably reflected in HRV [19]. So, in order to 
understand the properties of the fast noise, one would ideally remove the respiratory 
perturbation and consider the residual HRV which would then reflect fluctuations of the 
intrinsic dynamics of the heart control system. 

Consideration of the intrinsic activity of the heart control system on short-time 
scales is important for general understanding of the function of the cardio-vascular 
system, leads potentially to diagnostics of causes of arrhythmia involving problems with 
neuronal control [20], and can be a benchmark for modeling HRV. In this paper we 
present the results of an experimental study of the intrinsic dynamics of the heart 
regulatory system and discuss these results in the context of modelling the fast noise 
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Figure 1. i?i?-intervals for (a) normal (spontaneous) and (b) intermittent respirations. 
Respiration signals (arbitrary units) are shown by dashed lines. 
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Figure 2. (a) An ECG signal and (b) the corresponding HRV (RR intervals) signal. 
In (a) the R-peaks are marked by O ; the ECG signal is shown in arbitrary units. 



component. A number of open problems are identified. 
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2. Experimental results 

We analyse the dynamics of the control system in the absence of explicit perturbations 
by temporarily removing the continuing perturbations caused by respiration [figure QJb)]. 
To do so, we perform experiments involving modest breath-holding (apncea) intervals. 
Note that during long breath-holding the normal state of the cardiovascular system is 
significantly modified |21| . The idea of the experiments came from the observation that 
spontaneous apncea occurs during repose. Apncea intervals of up to 30 sec were used, 
enabling us to avoid either anoxia or hyper- ventilation [2"T] . 

Respiration-free intervals were produced by intermittent respiration, involving an 
alternation between several normal (non-deep) breaths and then a breath-hold following 
the last expiration, as indicated by the dashed line in figure H(b). The respiratory 
amplitude was kept close to normal to avoid hyper- ventilation, and there were relatively 
long intervals of apncea when the heart dynamics was not perturbed by respiration. It 
is precisely these intervals that are our main object of analysis. The durations of both 
respiration and apncea intervals were fixed at 30 sec. 

Measurements were carried out for 5 relaxed supine subjects, and they were 
approved by Research Ethics Committee of Lancaster University. Note that the 
measurements presented have been selected from a larger number of measurements to 
form a set recorded under almost identical conditions of time and duration, with the 
subjects avoiding either coffee or a meal for at least 2 hours beforehand. They were 4 
males and 1 female, aged in the range 29-36 years, non-smokers, without any history 
of heart disease. We stress that the aim of the current investigation was exploratory: 
to study typical behaviour of the internal regulatory system; we have not performed a 
large-scale trial of the kind widely used in medicine when a large number of subjects 
is necessary because of the need for subsequent statistical analysis of the data. The 
electrocardiogram (ECG) and respiration signals were recorded [T7j over 45-60 minutes. 
The ECG signals were transformed to HRV by using the marked events method for 
extraction of the RR- intervals which are shown in figure [2j 

Figure [1] shows i?i?-intervals found for the different types of respiration. It is 
evident that respiration changes the heart rhythm very significantly. Immediately 
after exhalation (b), there is an apncea interval where the heart rhythm fluctuates 
around some level. These fluctuations correspond to the intrinsic dynamics of the heart 
control system. It is clear from (a) that heart rate is continuously perturbed during 
normal respiration, whereas in (b) one can distinguish an interval of intrinsic dynamics 
corresponding to apncea. Thus, the jth interval of apncea is characterized by the time 
series {RRi}; here i — 1, 2 . . . labels the ith i?i?-interval. Finally, we form a set {RRiY 
for analyses by considering the set as realizations of a random walk and analyzing their 
dynamical properties as such. 

To reveal dynamics additional to -RR-intervals, the differential increments ARRi = 
RRi+i — RRi were analyzed. The differences between i?i?-intervals and their increments 
are illustrated in figure [3j Each apncea time-series {RRiY exhibits a trend that is 
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Figure 3. (a) i?i?-intervals and (b) increments ARR corresponding to apnoea intervals 
are shown. For convenience of presentation, the difference between a given value and 
the first value of each jth apncea interval is drawn in each case: RR\ = RRl — RR{ 
and ARR\ = ARR{ - ARR{. 



describable by the slope a of a linear function RR\ oc aj i, where % is a heart beat number 
and j marks jth apncea interval. The trend can be characterized by the distribution of 
slopes P(a) shown in figure H] (a). For all measurements the distributions are broad and 
their mean values differ from zero. Thus the non- stationary nature of HRV on short 
time-scales is clearly apparent. Note, that the distributions p(a) for the increments 
ARR are significantly narrower [figure H] (b)] and that they are very well fitted to a 
normal distribution; however, the mean values of the slopes differ from zero. 

Because the dynamics of i?i?-intervals is evidently non-stationary, we have applied 
detrended fluctuation analysis (DFA) [10] for estimation of the scaling exponents (3 
for the apncea sets {RRiY . In doing so, we adapted the DFA method [10] for short 
time-series and used non-overlapped windows (see Appendix for details). Because the 
time-series were short, time windows of length 4-15 i?i?-intervals were used to calculate 
(3. For all measured subjects, this procedure yielded values of (3 lying within the range 
(3 G (1.3 : 1.7), with a mean value of 1.45. If i?i?-intervals in the sets {RRiY are 
replaced by realizations of Brown noise (the integral of white noise) keeping same lengths 
of apncea intervals, then the calculation gives (3 = 1.46 ±0.07. Additionally, a surrogate 
analysis was performed for each subject by random shuffling of the time indices i of 
i?i?j-intervals, to confirm the importance of time-ordering of the i?i?-intervals. For 
each realization (set {RRiY), 100 surrogate sets were generated, 100 values of (3 were 
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Figure 4. Distributions of trend slopes P(a) of the sets (a) {RRiY and (b) {ARRi}^ 



obtained, and the mean value j3 s was calculated. Values of (3 S for the surrogate sets 
lie in the same limits as those for the original sets, but with a small bias between (3, 
calculated using original sets, and f3 s (see the Appendix for values of (3 and (3 S ). It means 
that one can see a correlation between i?i?-intervals, but that it is weak. Summarizing 
the DFA results, we can claim that the scaling exponent (3 is similar to that for free 
diffusion of a Brownian particle, but there is nonetheless some correlation between the 
i?i?-intervals. We also applied aggregation analysis [19] in a similar manner and arrived 
at qualitatively the same conclusion. Note that in the contrast to the initial idea of the 
DFA and aggregation analyses, which were used for revealing long-range correlations 
in time series, we have used these approaches to analyse the diffusion velocity because 
they can cope with trends. Long-range correlations cannot be revealed in the described 
measurements. 

To estimate the strength of the correlation, stationary time-series of the increments 
{ARRiY were considered. The autocorrelation function p(k) was calculated 
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Here RR ri = ARR? — (ARRi); the brackets () denote calculation of the mean value; % and 
j correspond to the heart beat number and apncea interval respectively, k = 0, 1, . . ., 
m? is the number of increments ARR in the j th apncea; N is the total number of 
apncea intervals. Figure [5] presents examples of autocorrelation functions. One of 
them has pronounced oscillations. An approximation of p(k) by the function p a (k) = 
exp(— jk) cos(27rfiA;) demonstrates that oscillations occur with frequency near 0.1 Hz, 
presumably corresponding to myogenic processes [TF] or (perhaps equivalently) to the 
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Figure 5. Examples of the autocorrelation function p(k) (a) with and (b) without 
an oscillatory component. The crosses indicate p(k) calculated on the basis of the 
increments ARR. The solid line corresponds to the approximating curve p a (k) — 
cxp(— 7&) cos(2irflk). 



Mayer wave associated with blood pressure feedback [22j [23] . Further investigations via 
the parametrical spectral analysis for each apncea interval show that these oscillations 
are of an on-off nature, i.e. observed for parts of the apncea intervals, and not in all 
of the measurements as can be seen in figure [5] (b). Examples of apncea intervals 
with and without oscillations are shown in figure [6j When an oscillatory component is 
present then its contribution to p(k) is much weaker than the contribution of the noisy 
component. The latter is characterized by a very short memory as demonstrated by fast 
decay of p(k). 

The properties of ARR can also be characterized by the probability density function 
P(ARR) shown in figured (a). Figure (b) shows the probability density function P(RR) 
of RR- intervals for comparison. Following [24] , the a-stable distribution has been widely 
used to fit the distribution of increments ARR, and strongly non-Gaussian distributions 
were observed |24j . We perform a similar fitting applying special software [25]. Since 
the distributions P(ARR) are almost symmetrical, our attention was concentrated on 
the tails, which were characterized by a stability index a G (0,2]. The case of a = 2 
corresponds to a Gaussian and, if a < 2, the tails are wider than Gaussian. Fitting to 
our results yields a stability index a G (1.8 : 2), and the goodness-of-fit test (modified 
KS-test taking into account the weight to the tails [25]) supports the fitting. Note that, 
although the autocorrelation function p(k) cannot be used for the theoretical description 
of an a-stable process [26] . p(k) is nonetheless applicable for finite time-series. 

If we consider the same length of realization using a Gaussian random variable 
instead, we find a = 1.99 ± 0.01. It means that the calculations of a are very robust. 
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Figure 6. Examples of apnoea intervals with (a) and without (b) oscillation of HRV. 
The circles correspond to the values of the increments ARRi and the solid lines 
connecting points are guides to the eye. The dashed lines in figure (a) are added 
to reveal oscillations. The middle and upper ARRi time-series are shifted by 0.1 and 
0.2 (sec) accordingly. 

In addition we carried out a stability test and it too supported the fitting results. The 
obtained values of a G (1.8 : 2) differ significantly from the previously reported values 
a G (1.5 : 1.7) for 24h time-series of i?i?-intervals [24J. 

Combining all the results, we conclude that the short-time dynamics of RR- 
intervals can be described as a stochastic process with stationary increments. This 
type of stochastic processes was discussed by A. N. Kolmogorov [27] and applied to the 
description of a number of different problems (see e.g. [281 1291 130] for further details). 
So, HRV during apncea interval cab be presented in the following form 

RRi = RRi—i + /ARRi, (2) 

where ARRi is a stationary discrete time stochastic process. Note that the DFA 
calculation excludes a linear trend, which is taken into account in Eq. (j2j) as non- 
zero mean value of the increments, fij = (ARRi)j; in general case, \ij is a random 
function of jth apncea interval. If one represents i?i?-intervals as a sum of the linear 
trend and a random component: 

RRi = fij i + (3) 
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Figure 7. (Color online) Normalized probability density functions (a) P(ARR) of 
increments of i?i?-intcrvals and (b) P(RR) of i?i?-intervals. In (a) the full (blue) and 
dashed (red) lines are Gaussian and stable distributions, respectively, fitted to the 
data. The insets show the same distributions plotted with logarithmic ordinate scales; 
the circles correspond to P(ARR). The stable distribution in (a) is characterized by 
a = 1.86. In (b) the full (blue) line is a Gaussian distribution fitted to the data. 

then £j corresponds to the non- stationary process (J2J) with zero mean value of increments. 
In other words, the superimposed random component of HRV during apncea intervals 
is described by a non- stationary random process. 

Increments ARRi are characterized by a random a-stable process of short memory, 
with a weak intermittent oscillatory component of frequency ~0.1 Hz. In the zeroth 
approximation the increments can safely be represented by an uncorrelated Gaussian 
random process but, in the next approximation, a weak correlation must be included, 
allowing for an intermittent oscillatory component, and for weak non-Gaussianity of 
the distribution of increments ARR. These additions reveal, on the one hand, that 
the previously reported observation of a non-Gaussian distribution of increments [21] 
is a property of the intrinsic heart rate regulatory system, but on the other hand, 
that the scaling ranges of the stability index a differ significantly in the presence 
or absence of external perturbations (including respiration) acting on the regulatory 
system. Consequently an explanation of the scalings reported in [2U [10] should include 
analyses of the effect of external perturbations and respiration, and not an analysis of 
heart rate alone. 

3. Discussion 

3.1. Non-stationarity of RR-intervals during apnoea 

The results presented indicate that there is no firm set point for the heart control system, 
and that the heart rhythm exhibits diffusive behaviour. The slowest dynamics can be 
described by a linear trend during apncea intervals and its presence can be treated as 
a slow regulatory/adaptation component of the control system. The presence of the 
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slow time-scales is an established property of HRV [21] and their presence, even in the 
absence of the respiratory perturbation, can be interpreted as an expected property. 

On short time-scales of order several seconds, HRV shows a diffusive dynamics too. 
It can be interpreted in two ways. One possibility is that the control system does not 
firmly trace the base (slow) rhythm, because in case of tracing, short time-fluctuations 
should "jump" around the base rhythm and, consequently, be stationary. Such a picture 
corresponds to zero action of the control system if the heart rate is in a "safe" (for the 
whole cardiovascular system) interval, e.g. RR G [RRi ow : RRhigh}- Another possible 
explanation could be that the control system is tracing the base rhythm but the short- 
time fluctuations have a non-stationary character. It is natural to expect that there 
could be other possible explanations, and additional investigations are needed to reach 
an understanding of the diffusive dynamics on short-time scales. 

In section 2 it is suggested that we should consider the non-stationarity and diffusive 
dynamics of i?i?-intervals within the framework of a stochastic process with independent 
increments. It allows one to consider i?i?-intervals as realizations of the so-called auto- 
regressive process that is widely used in time-series analysis [32]. It means that the 
direct spectral estimation of i?i?-intervals, currently used as one of the basic techniques 
[31], is not applicable here and that one must use the theory of stochastic processes 
with stationary increments for their spectral decomposition [30]. If in the presence 
of respiration, the short-time stochastic component of HRV preserves non-stationarity 
then spectral estimation based on i?i?-intervals is not correct, and increments must be 
used instead. Note, that the properties of short-time fluctuations in the presence of 
respiration are far from being completely understood. 

3.2. Non-Gaussianity and correlations of increments ARR 

The theories of both stochastic processes with stationary increments and of auto- 
regressive analysis place some limitations on the analysed time-series. The first approach 
requires the existence of finite second-order momenta, whereas the second approach 
assumes uncorrelated statistics of increments. Formally, however, non-Gaussianity of the 
increments distribution means that the second-order momenta do not exist [26], but non- 
Gaussianity can still be incorporated into the auto- regressive description [33]. And vice 
versa, the presence of correlations in the increments dynamics requires a modification of 
the standard auto-regressive approach, and it is one that can be incorporated naturally 
into the general theory. In the current investigation we ignore these issues. We calculate 
the auto-correlation function and use model (j2J), because the finite length of the time- 
series guarantees the existence of the second-order momenta, and the simplicity of ([2]) 
means that the inclusion of the correlations is a trivial extension. 

Our consideration has the formal character of time-series analysis because we do not 
incorporate any preliminary information about the possible dynamics of i?i?-intervals. 
The analysis is based on the use of a set of relatively short time-series, a fact that defines 
our choice of simple statistical measures. One cannot exclude the possibility that the use 
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of other approaches to such data might provide additional insight into HRV dynamics. 
For example, the fractional Brownian motion approach [191 EI] and the theory of discrete 
non- stationary non-Markov random process [35] represent different paradigms, which are 
based on assumptions about the origin of the data. Note, that despite a long history of 
developing the approaches and their applications, the approaches of fractional Brownian 
motion and of stable random process are not standardized tools, whereas the approach of 
non-Markov random process is not so popular. There is no definite recipe for choosing 
a set of measures which can uniquely specify (or provide a good description of) the 
properties of a renewal (discrete time) stochastic process. 

3.3. Modelling 

Another way of attempting to understand the results is to try to reproduce the observed 
data properties from an appropriate model. In the context of our experiments, the 
modelling should consist of a simulation of the electrical activity of sinoatrial node 
(SAN) where the heart beats are initiated. For modelling, one option is to use a bottom- 
up approach, which is currently a very popular technique within the framework of the 
complexity paradigm. In fact, available SAN cellular models allow one to incorporate 
many details of physiological processes like the openings and closures of specific ion 
channels [36]. However, despite the complexity of the models (40-100 variables) many 
important features are still missed. For example, the fundamentally stochastic dynamics 
of ion channels is represented by equations that are deterministic. Heterogeneity of the 
SAN cellular locations and intercell communications are among other important open 
issues [371 EE] • 

An alternative option is the top-down approach using integrative phenomenological 
models. In contrast to detailed cell models, a toy model of the heart as a whole unit can 
be developed. It is known that an isolated heart, and a heart in the case of a brain-dead 
patient |39j demonstrate nearly periodic behaviour. So, it is reasonable to assume that 
the observed HRV is induced by the neuronal heart control system, which is a part of 
the central nervous system. The control system includes a primary site for regulation 
located in the medulla [40J , consisting of a set of neural networks with connections to the 
hypothalamus and the cortex. The control is realized via two branches of the nervous 
system: the parasympathetic (vagal) and the sympathetic branches. Although many 
details of the control system are still missing [HI [20], it is currently accepted that the 
vagal branch operates on faster time scales than the sympathetic one, and that each 
branch has a specific co-operative action on the heart rate and the dynamics of SAN 
cells. 

Let us consider an integrate-and-fire (IF) model as a model of a SAN cell in the 
leading pacemaker. These cells are responsible for initiating the activity of SAN cells 
and, consequently, that of the whole heart [38] • The dynamics of the IF model describes 
the membrane potential U(t) of the cell by the following equations 



dU 
~dt 



T 



1 



if U(t)<U t 



(4) 
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U(t) = U r , t* = t if U(t) = U t and ^ > 0, (5) 

Here 1/r defines a slope of integration, Ut is the threshold potential, U r is the resting 
(hyperpolarization) potential; the time t* corresponds to the cell firing, and it is the 
difference between two successive firings that determines the instantaneous heart period 
or i?i?-interval, RRi = t* — t*_ 1 . It is known [ID] that increasing sympathetic activity 
with a combination of decreasing vagal activity leads to an increase in the heart period, 
and vice versa. Direct stimulation of the sympathetic branch leads to an increase of 
the integration slope 1/r and a lowering of the threshold potential Ut, whereas vagal 
activation has the opposite effects, and additionally, lowers the resting potential U r . 
Thus, the neuronal activities can be taken into account as modulations of the parameters 
of the IF model (jlj). For reproducing HRV during apncea, therefore, it is enough to 
present any of the parameters r, U t or U r as a stochastic variable of the form (j2J), 
for example, U t (t*) = Utit*^) + £ i( where & are random numbers having the stable 
distribution. 

However, the use of more realistic (than IF) models with oscillatory dynamics, for 
example Fitzhugh-Nagumo [12] or Morris-Lecar [13] models, makes the reproduction of 
the experimental results a more difficult but interesting task. Currently it is unclear 
whether it is possible to obtain a stable distribution of increments by consideration 
of the Gaussian type of fluctuations alone, or whether one should use fluctuations 
characterizing by a stable distribution. This point demands further investigation. 

4. Conclusion 

In summary, our experimental modification of the respiration process reveals that the 
intrinsic dynamics of the heart rate regulatory system exhibits stochastic features and 
can be viewed as the integrated action of many weakly interacting components. Even on 
a short time scale (less then half a minute) the heart rate is non-stationary and exhibits 
diffusive dynamics with superimposed intermittent ~0.1 Hz oscillations. The intrinsic 
dynamics can be described as a stochastic process with independent increments and 
can be understood within the framework of many-body dynamics as used in statistical 
physics. The large number of independent regulatory perturbations produce a noisy 
regulatory background, so that the dynamics of the regulatory rhythm is close to classical 
Brownian motion. However there are indications of non-Gaussianity of increments 
and weak but important correlations on short time-scales. The reproduction of these 
features, especially the non-Gaussianity property, is an open problem even in simple toy 
models. 

These results are important both for understanding the general principles 
of regulation in biological systems, and for modeling cardiovascular dynamics. 
Furthermore, the results presented may possibly lead to a new clinical classification 
of states of the cardiovascular system by analysing the intrinsic dynamics of the heart 
control system as suggested in [20j. 
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Table Al. Data for each subject. (RR) is the mean heart rate during apncea. J 
is the total number of apnoea intervals. /3 is the DFA scaling exponent calculated for 
the apncea set {itRjp. (3 S is the mean value of the DFA scaling exponent calculated 
for surrogate data, which were generated by random shuffling of the time indices i 
of i?i?i-intervals. b is the scaling exponent of the aggregation analysis. 7 and f2 
correspond to the values of parameters for the function p a (k) = cxp(— 7&) cos(27rf2fc) 
which approximates the autocorrelation function p(k). a is the stability index of the 
distribution P(ARR). 
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SI 


1.01 


45 


1.39 


1.47 


1.86 


0.81 


0.17 


1.83 


S2 


0.77 


4G 


1.46 


1.44 


1.83 


0.21 


0.09 


1.95 


S3 


1.10 


47 


1.43 


1.53 


1.96 


1.01 


0.22 


1.79 


S4 


0.75 


47 


1.58 


1.60 


1.91 


0.15 


0.08 


1.90 


S5 


0.91 


GO 


1.42 


1.48 


1.82 


0.28 


0.13 


1.86 
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Appendix 

Some details of the measurements and calculations are summarized in this section. 

The ECG was measured by standard limb (Einthoven) leads and the respiration 
signal was measured by a thoracic strain gauge transducer. The signals were digitized 
by a 16-bit analog-to-digital converter with a sampling rate of 2 kHz. The ECG and 
respiration signals were recorded over 45-60 minutes and time locations of .R-peaks in 
the ECG signals were defined and time intervals between two subsequent -R-peaks (the 
so called i?i?-intervals) are used to form HRV signal. 

Respiration-free intervals were produced by the intermittent respiration, involving 
an alternation between normal breaths and apncea intervals. The durations of both 
normal breaths and apncea intervals were fixed at 30 sec. The respiration signal 
was used to identify apncea intervals. Finally, the set of time-series of RR-interval 
{RRiY was formed for each subject; here i — 1,2 . . . labels the ith i?i?-intervals, and 
j = 1,2. . . labels the jth interval of apncea. For each interval of apnoea, time series 
of the differential increments ARRi = RRi+\ — RRi were produced and they also form 
a set {ARRiY for each subject. The number of i?i?-intervals in each apncea interval 
is different, depending on the heart rate of the subject. The total number of apncea 
intervals also differ for each subject. The mean heart rate (RR) during apncea intervals 
and the total number J of intervals for each measured subject are presented in table 
EH 

For the application of the DFA and aggregation analyses we adapted the approaches 
described in pi)] and [TH], respectively, to treat the available sets of short time series 
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{RRiY. 

The DFA exponent (3 was calculated in the following way. First, the initial set 
{RRiY was transformed to another set {y(k)Y by the following expression: 

k 

y(k) = J2RXii (A.l) 

8=1 

where k = 1 . . . Mj and Mj is the number of i?i?-intervals for jth apnoea interval. 
For each length n = 4, . . . 15 of time window a set of linear trends {y n {k)Y was 
calculated (see [JU] for details), where y n {k) = k ■ + 6^, m = 1,... [Mj/n\, 
\_x\ = max{n G 7L\n < x} is the floor function of x. Then a set of scaling function 
{F(n)Y was calculated for each value of n by use of the expression 

F{n) = Y J W)-y n {k)}\ (A.2) 

k 

where Nj = \_Mj/n\ ■ n. Further the scaling functions F{n) were calculated as 



F(n) 



\ 



3=1 



where J is the number of apnoea intervals for the given subject, iV = $^/=i Nj. Finally, 
the scaling exponent (3 was determined as a slope of the function log[F(n)] oc /31og(n) 
(see figure IAT1 (a)). The values of (3 for the different subjects are shown in table [ATI 

The aggregation analysis consists of three steps and the final result is the scaling 
exponent b. The first step is the creation of a set of aggregated time series {z m {k)y for 
different m — 1, ... 10: 

k+m 

Zm{k) = J2 RR i> ( A ' 4 ) 

i=k 

where k — 1, . . . Mj — m. Then a realization z m (k) was formed from the set {z m (k)Y: 
Zm{k) = {z m (k)y = {zmik)} 1 , . . . {z m {k)Y ■ The second step includes the calculation 
of the mean value /i(m) and variance a(m) of the time-series z m (k): 

^ M M 
k=l k=l 

where M is the whole length of time series z m (k). The slope b of the function 
log[er(m)] oc 61og[/i(m)] was calculated in the third step (see figure IAT1 (b)). The values 
of b for each subject are shown in table IA1I 

To verify the robustness of the calculations of exponents (3 and b we have performed 
calculations with the same number of RR-intervals as well as the same structure of 
apnoea intervals but by using realizations of Brown noise generated by computer. In 
other words, in the procedures described above we replaced {RRiY by where 
Wi = Wi_i + 0.2 • £j for i = 2, • • • Mj, W\ = RRi, and & are random numbers having 
the normal distribution with mean zero value and unit variance; the numbers £j are 
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Figure Al. (a) The scaling function F(n) (circles) and its approximation (dashed 
line) by F(n) oc nP (/? — 1.39) are shown, (b) The dependence (circles) of the variance 
cr(m) on the mean value fi(m) for to = 1, ... 10 and its approximation (dashed line) 
by aim) oc n(m) b (b = 1.82) are shown. 



different for different jth intervals of apneea. We performed 100 calculation of (3 and b 
for different sets {H^p for each subject. Theoretical values of (3 and b for the Brown 
noise are /? = 1.5 and 6 = 2 correspondingly. The calculations with Brown noise gave 
(3 = 1.46 ± 0.07 and b = 1.84 ± 0.04. Here data were merged for all subjects and are 
presented in the form of a mean value ± its standard deviation. It means that there is a 
systematic error related to the length and data structure, a general error of calculation 
in respect to the theoretical values for (3 is 1.5 and for b is 2.0. However the standard 
deviations of the calculated values are rather small and, consequently, we can conclude 
that our calculations of (3 and b are robust. 
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